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sJ ■ Abstract 

B\ 

The ratio of the photon mediated production or annihilation rates of spin triplet 
and spin singlet heavy quarkonium states is computed to the next-to-next-to- leading 
logarithmic accuracy within the nonrelativistic renormalization group approach. 
The result is presented in analytical form and applied to the phenomenology of tt, 
bb and cc systems. The use of the nonrelativistic renormalization group consider- 
ably improves the behaviour of the perturbative expansion and is crucial for accu- 
rate theoretical analysis. For bottomonium decays we predict r(%(15) — > 77) = 
0.659 ± 0.089(th.)JIo:oi8(<fas) ± 0.015(exp.) keV. Our results question the accuracy 
of the existing extractions of the strong coupling constant from the bottomonium 
annihilation. As a by-product we obtain a novel result for the ratio of the ortho- 
and parapositronium decay rates: the corrections of order a 4 In a and a 5 In a. 

PACS numbers: 12.38.Cy, 13.20.Gd 



1 Introduction 

Annihilation of a heavy quarkonium into leptons, photons or light hadrons, as well as 
its production in e + e~ or 77 collisions, has been a subject of numerous investigations 
starting from earliest applications of perturbative quantum chromodynamics (QCD) pQ 
and has by now become a classical problem. The annihilation of the nonrelativistic bound 
state possesses a highly sophisticated multi-scale dynamics and, not surprisingly, was the 
original testing ground for the effective field theory of nonrelativistic QCD (NRQCD) 
[2f3] . On the phenomenological side, the ratio of the bottomonium annihilation rates into 
leptons and light hadrons is currently used to determine the value of the strong coupling 
constant a s [4J. Furthermore, the study of ti threshold production at a future Linear 
Collider could even allow us to probe Higgs-boson-induced effects [5]. The calculation of 
the high-order QCD corrections is mandatory to provide accurate theoretical predictions 
for the production and annihilation rates, however, this is a challenging theoretical prob- 
lem. For most of the channels, the perturbative corrections are only known to one-loop 
accuracy [6]. The complete 0(a 2 s ) correction is known only for the one- and two-photon 
mediated processes [7< 8ll9lll0|lll|ll2lll3|ll4] . The perturbative series, however, shows a slow 
convergence and the results suffer from a rather strong dependence on the normalization 
scale v of ct s [y\ The resummation of the large logarithms of the heavy quark velocity to 
all orders in a s has been advocated as a tool to improve the behaviour of the perturba- 
tive expansion for tt threshold production [15]. Currently, the complete next-to- leading 
logarithmic (NLL) approximation for the production and annihilation rates is available 
[T6p7] . The first attempt to go beyond the NLL approximation [15] suggested a very 
good convergence of the logarithmic expansion. In particular, an accuracy of 2-3% was 
claimed for the cross section of ti threshold production. However, subsequent calculations 
of some next-to- next-to-leading logarithmic (NNLL) terms |18j . which had not been taken 
into account in Ref. [15] , casted serious doubts on this optimistic estimate. Thus, the full 
calculation of the NNLL corrections, which still remains elusive, is unavoidable to draw 
definite conclusions. 

In this paper we derive the complete NNLL result for the spin dependent part of 
the heavy quarkonium production annihilation rates which includes the terms of the 
form a™ +2 ln n a s for all n. The structure of the paper is as follows. In the next section 
we recall the basic ingredients of the effective theory description of the production and 
annihilation rates. In Sect. [3] we derive the nonrelativistic renormalization group (NRG) 
equation and in Sect. H] we present its solution. The result is applied to ti threshold 
production, bottomonium and charmonium phenomenology in Sect. EJ Sect. M contains 
our conclusions. Explicit results for the potentials and the analytical solution of the 
equations are given in Appendices A and B, respectively. Appendix C includes the new 
results for the high-order logarithmic corrections to the positronium decay rates. 



2 Effective theory description of the production and 
annihilation rates 

The perturbative dynamics of the nonrelativistic bound state is characterised by three well 
separated scales: the heavy quark mass m q (the hard scale), the bound state momentum 
vm q (the soft scale), and the bound state energy v 2 m q (the ultrasoft scale), where v ex 
a s <C 1 is the velocity of the nonrelativistic heavy quark for the approximately Coulombic 
bound state. A systematic evaluation of the corrections to the bound state parameters 
is based on the effective field theory concept of the scale separation [2], which is at the 
heart of the recent progress in the perturbative QCD bound-state calculations. 

The one-photon mediated processes are induced by the electromagnetic current j^. 
Its space components have the following decomposition in terms of operators constructed 
from the nonrelativistic quark and antiquark two-component Pauli spinors ip and x [3] : 

j = c v (u)^ax + d -^^aD\ + . . . , (1) 

where v is the renormalization scale, D is the covariant derivative, a is the Pauli matrix, 
and the ellipsis stands for operators of higher mass dimension. The Wilson coefficients c v 
and d v represent the contributions from the hard modes and may be evaluated as a series 
in a s in full QCD for free on-shell on-threshold external (anti) quark fields. We define it 
through 

cM = £(^)>(,), 4°) = !, (2) 

and similarly for other coefficients. 

The operator responsible for the two-photon S"-wave processes in the nonrelativistic 
limit is generated by the expansion of the product of two electromagnetic currents and 
has the following representation [3] 

77 = c 7 »^ x + d -^^D\ + . . . , (3) 

which reduces to the pseudoscalar current in the nonrelativistic limit. The one-loop 
corrections to the hard coefficients are known for quite a long time [19.20J 
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t" = -2C, 



where Cp = (N% — 1)/(2N C ), N c = 3. Let us define the spin ratio for the production and 
annihilation of heavy quarkonium Q as 



a(e+e" -> Qj^S^) _ r(Q(n 3 g x ) -> e + e 
<r( 7 7 - QCn^o)) " T{Q{riS Q ) - 77 ) 



n i = 1777 . T^TZTB ^ = WT^TZTB \ . ...a • ( 5 ) 



The effective theory expression for the spin ratio reads 
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where Q q is the quark electric charge, c s (v) = c„(z/)/c 77 (z/), 4 } n' p \ r ) are the spin triplet 
(vector) and spin singlet (pseudoscalar) quarkonium wave functions of the principal quan- 
tum number n. The wave functions describe the dynamics of the nonrelativistic bound 
state and can be computed within potential NRQCD (pNRQCD) [21]. The latter is the 
Schrodinger-like effective theory of potential (anti) quarks whose energy scales like m q v 2 
and three-momentum scales like m q v , and their multipole interaction to the ultrasoft glu- 
ons [22|23] . The contributions of hard and soft modes in pNRQCD are represented by 
the perturbative and relativistic correction to the effective Hamiltonian, which is system- 
atically evaluated order by order in a s and v around the leading order (LO) Coulomb 

^£(0) 2 = C 3 F a 3 s m 3 /(8rrn 3 ). Let 
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The LO Coulomb wave function reads 
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The next-to-leading (NLO) contribution pW vanishes since the corrections to the wave 
functions are spin-independent at this order. 

Starting from 0(af), the hard coefficients are infrared (IR) divergent. This spurious 
divergence arises in the process of scale separation and is canceled against the ultraviolet 
(UV) one of the effective-theory result for the wave function at the origin. A powerful 
approach to deal with such divergences has been developed in Refs. [2"4l25f26lT2"7] . It 
is based on dimensional regularization and the interpretation of the formal expressions 
derived from the Feynman rules of the effective theory in the sense of the threshold 
expansion [28.29J. This provides a factorization of the contributions from different scales. 
In the MS the two-loop corrections to c v are known in analytical form and given by [HI 
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where (3q = 11Ca/3 — 4n/Tp/3, Ca = N c , T F = 1/2, rii is the number of light-quark 
flavours, £(3) = 1.202057. . . is the value of Riemann's £ function, and a s is renormalized 
in the MS scheme. For the two-photon processes the two-loop correction is known in 
semi-numerical form Hi 
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where the (7(7 — > 77 contribution induced by a light-fermion box was estimated to be small 
and was not included in the result. Note that the above result depends on the definition 
of the nonrelativistic axial current. The next-to-next-to-leading (NNLO) correction to 
the wave functions ratio in MS scheme reads 1141301 
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where ^f n (x) = d n \nT(x)/dx n , V(z) is Euler's T function, and 7^ = 0.577216 ... is Euler's 
constant. However, for the renormalization group analysis it is more convenient to use the 
hard matching scheme where the nonlogarithmic part of the divergent two-loop potential- 
potential contribution, which appears in the time-independent perturbation theory for the 
nonrelativistic wave function, is shifted to the hard coefficient. The NNLO corrections in 
this scheme are obtained from the MS result by the following shift 
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3 Renormalization group evolution of c s 

In general, the logarithmic corrections originate from logarithmic integrals over virtual 
momenta ranging between the scales and reveal themselves as the singularities of the 
effective theory couplings. The renormalization of these singularities allows one to derive 
the equations of the nonrelativistic renormalization group (NRG), which describe the 
running of the effective theory couplings (Wilson coefficients), i.e. their dependence on 
the effective theory cutoffs. The solution of these equations sums up the logarithms of the 
scale ratios. To derive NRG equations necessary for the NNLL analysis of the decay rates 
we rely on the method developed in Ref. [16] where, in particular, the correct NLL result 
for the decay rates has been obtained for the first time (see also Ref. [IT]). The NRG 
equations express the dependence of the effective theory coupling constants on the IR 
cutoff. However, they are derived by studying the UV divergences of the effective theory 
perturbative expressions. In general, one has to consider the soft, potential and ultrasoft 
running of the effective theory coupling constants. We denote the corresponding cutoffs 
as v s , Up and v us , respectively. The soft running is associated with the divergences of the 
NRQCD perturbation theory for the potential while the potential running corresponds to 
the divergences of the time-independent perturbation theory for the nonrelativistic Green 
function in pNRQCD. As it was first realized in Ref. [31], v us and v v are correlated and 
the relation between them can be given by v us = vi/m q . The matching to the hard 
contribution is performed at a generic scale Uh ~ m q . 

To NNLL approximation the hard coefficient c s itself has only potential running. We 
compute the corresponding anomalous dimensions by inspecting the UV singular be- 
haviour of the three-loop effective theory diagrams computed in Refs. [32,30] (see also 



[33|34f35| ). The NRG equation in the hard matching scheme is found to be 
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Here av s ( u ) — a s (u) + . . . is the strong coupling constant defined through the perturbative 
potential between the static quark and antiquark in colour-singlet state 

2 47rC F a Vs (k 2 ) 
v c{k*) = p , (13) 

and D s i (y) = a s {v) + . . . is the Wilson coefficient of the spin-flip potential in pNRQCD 
Hamiltonian 
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with S = {(T\ + cr 2 )/2 being the total spin operator. In Eq. ( Tl2|) we combine the soft and 
potential running by setting v s = v v = u, which is consistent to the order of interest. Note 
that the scheme used in this paper differs from the standard one used for the evaluation 
of the Feynman integrals in dimensional regularization in Ref . [30] . 

Let us now discuss the structure of the NRG equation (TT2|) in more detail. The 
first term on the right-hand side originates from the diagram (a) in Fig. [TJ The leading 
logarithmic (LL) soft running of ay a and D^J is responsible for the NLL evolution of 
c s . To get the NNLL result for the hard coefficient we need the NLL running of this 
quantities. The one of ay s is well known and can be found e.g. in Ref. [IB]. For Daa it 
has been recently obtained in Refs. [36|37] . The analytical result for the NLL evolution 
of the spin-flip potential is presented in Ref. [37] where, however, some details of the 
analysis have been skipped. This gap is filled in Appendix A, where explicit results for 
the contributing potentials are presented. 

The second term on the right-hand side of Eq. (TT2"|) originates from the diagram (b) 
of Fig. [TJ It starts to contribute to c s at NNLL. Thus the LL expression for DL \ {v) can 
be used, which is given by [T6|38] 

^_ D g) M = f?M 7 w c | H , 7 w = _| + 7 Ca . (15) 

ahiv ' 7i 2 4 

cf is the effective Fermi coupling, which is needed with LL accuracy, and reads c F {i>h) — 
z~ Ca , where z = (a s (v)/a s (v h )) l/Po . 

The initial condition for the NNLL solution of the NRG evolution is fixed by the 
two-loop value of c s (uh) at a hard matching scale vu ~ "V The subsequent evolution of 
c s {y) down to v ~ a s m q resums the logarithms of the coupling constant. The form of the 
NRG equation and the matching conditions is scheme dependent. A change of the scheme 
Csi^h) —* c fK u h) = cf\u h ) — 5c s tt 2 C f requires an additional full derivative term to be 



(a) (b) 

Figure 1: The pNRQCD diagrams contributing to the NNLL running of c s . The shaded oval 
represents the soft running of the spin-flip or Coulomb potential and j is nonrelativistic vector 
or pseudoscalar current. 



added to the right hand side of Eq. (Jl2j) 



(16) 
where Sc s specifies the scheme and the ellipses indicated higher orders in a s . In particular 
we have [5c s ] MS = —7/9. Note that the coupling cp appears in Eq. (TT2|) in combination 
with the factor 7W. 

To get the NNLL approximation for the production and annihilation rates we formally 
have to take into account also the LL running of d v and <i 77 . These hard coefficients, 
however, have identical LL ultrasoft running which cancels out in the spin ratio. 

4 Solution of the renormalization group equation 

The solution of the renormalization group equation (TT2l) is of the following form 

6.(1/) = c s(l/t)e « s WCM^WrrM^ . (17) 

The expression of T^ LL is known [TBIfPT] 

T NLL M = 27rC %(2A) - 7C A ) , Pn -2C A \ 2nC 2 F C A p 

Cs {) 3((3 -2C A y V 1 ) p M -2C A ) [n[Z ^ ^ 

The result for the NNLL function can also be obtained analyticalljo and can be cast in 
the form 

19 

rr>) = ^ 2 E^W, (is) 



1 For the practical calculation we use v us = v 2 jvh, which is sufficient for the accuracy of our calculation. 



where the coefficients A4 and functions fi(z) are given in Appendix B. To get the NNLL 
approximation for the spin ratio one has to take into account the NNLL contribution to 
the wave function, which is given by 

pr» = 1 + ^ (^>)4i(-)) LL pi 2) (-) • (20) 

It cancels the NLL //-dependence of c s . Eqs. (fTTj) and (1201) give the NNLL approximation 
to the spin ratio. By expanding the resummed expression we reproduce the known results 
for the 0(a 2 s ) [H], 0(a 3 s \n 2 a s ) [32], and 0(a 3 s \na s ) [313] terms. After including the 
one-photon annihilation contribution, the Abelian part of our result reproduces the spin- 
dependent corrections of order a 3 In a to the positronium decay rates [331135] . The use 
of the NRG allows us to derive the higher order logarithmic corrections for positronium. 
The explicit results for the 0(a 4 ln 2 a) and 0(a 5 ln 3 a) terms are given in Appendix C. 

5 Heavy quarkonium phenomenology 

For the numerical estimates, we take m c = 1.5 GeV, m^ = M (T) /2 and m t = 175 GeV, 
which is sufficient at the order of interest. Furthermore, we take a s [Mz) a s an input and 
run with four-loop accuracy down to the matching scale v^ to ensure the best precision. 
Below the matching scale, the running of a s is used according to the logarithmic precision 
of the calculation in order not to include the corrections beyond the NNLL accuracy. By 
the same reason we expand the decay rates ratio in a s (uh) up to the NNLL accuracy. In 
Figs. El [3] and HI the spin ratio is plotted as a function of v in the various logarithmic 
and fixed-order approximations for the (hypothetical) toponium, bottomonium and char- 
monium ground states, respectively. As we see, in the second order the convergence and 
stability of the result with respect to the scale variation is substantially improved if one 
switches from the fixed-order to the logarithmic expansion. We want to remark that the 
v dependence of the NLL approximation is slightly worse than at NLO. This is due to the 
artificially small v dependence at NLO which is likely due to the fact that at this order 
only the hard scale enters. 

Let us first consider the top quark case. The relatively large top quark width smears 
out the Coulomb-like poles below the threshold leaving a single well-pronounced resonance 
with the properties mainly determined by the "would be" toponium ground state param- 
eters. Thus our result can be applied to the ratio of the cross sections of the resonance 
e + e _ — > ti and 77 — > ti production. As one can see in Fig. [2] the logarithmic expansion 
shows perfect convergence and the NNLL correction vanishes at the scale v w 13 GeV, 
which is close to the physically motivated scale of the inverse Bohr radius a s m t /2. For 
illustration, at the scale of minimal sensitivity, v = 20.2 GeV, we have 

a^e+e-^tt) = 1 _ _ 

<T M (77->tt) 3 Qt 

Note that the perturbative expansion for the ground state energy, which is known to 
0(a 3 ) [39], shows as similar nice property. However, it is not clear if the nice behaviour 
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Figure 2: The spin ratio as the function of the renormalization scale v in LO=LL (dotted 
line), NLO (short-dashed line), NNLO (long-dashed line), NLL (dot-dashed line), and 
NNLL (solid line) approximation for the (would be) toponium ground state with u^ = m t . 
For the NNLL result the band reflects the errors due to a s {M z ) = 0.118 ± 0.003. Note 
that for the vertical axis the zero is suppressed. 



of the logarithmic expansion also holds for the spin-independent part of the threshold 
cross section. A possible problem is connected to the ultrasoft contribution, which is 
enhanced by the larger value of a s at the ultrasoft scale. Whereas it is suppressed in the 
spin ratio by the fifth power of a s , for the spin- independent part it already contributes at 
0(a\) and can destabilize the expansion. 

For bottomonium, the logarithmic expansion shows nice convergence and stability (c.f. 
Fig. [3]) despite the presence of ultrasoft contributions with a s normalized at a rather low 
scale v 2 jvfih- At the same time, the perturbative corrections are important and reduce 
the leading order result by approximately 41%. For illustration, at the scale of minimal 
sensitivity, v = 1.295 GeV, we have the following series: 



T(T(1S) -> e + e~; 
T( Vb (lS) -+ 77) 
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Figure 3: The spin ratio as the function of the renormalization scale v in LO=LL (dotted 
line), NLO (short-dashed line), NNLO (long-dashed line), NLL (dot-dashed line), and 
NNLL (solid line) approximation for the bottomonium ground state with v^ = m&. For 
the NNLL result the band reflects the errors due to a s (M z ) = 0.118 ± 0.003. 

In contrast, the fixed-order expansion blows up at the scale of the inverse Bohr radius. 

So far we have discussed the perturbative corrections to the Coulomb-like quarkonium. 
However, in contrast to the tt system, for bottomonium nonperturbative contributions can 
be important. In our case the interaction of the quark- ant iquark pair to the nonperturba- 
tive gluonic field is suppressed by v through the multipole expansion. In every order of the 
multipole expansion the nonperturbative contribution is given by the convolution of the 
quantum mechanical Green function with a nonlocal nonperturbative gluonic correlator. 
In general, we know little about the structure of these quantities. If v 2 m q 3> Aq C d, it 
can be investigated by the method of vacuum condensate expansion [mPT] . The result- 
ing series, however, is not expected to converge well in our case and suffers from large 
numerical uncertainties [121H3] . Nevertheless, it is possible to estimate the nonrelativistic 
suppression of the leading nonperturbative effect. For the ratio of the decay rates the 
nonperturbative contribution emerges via the 0(v A ) chromomagnetic dipole interaction 



10 



& 0.6 
0.4 



__.-.-■-■-" 



0.2 




0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 

v (GeV) 



Figure 4: The spin ratio as the function of the renormalization scale v in LO=LL (dotted 
line), NLO (short-dashed line), NNLO (long-dashed line), NLL (dot-dashed line), and 
NNLL (solid line) approximation for the charmonium ground state with v^ = m c . For the 
NNLL result the lower (yellow) band reflects the errors due to a s (Mz) = 0.118 ± 0.003. 
The upper (green) band represents the experimental error of the ratio [I] where the central 
value is given by the horizontal solid line. 

and the interference between the 0(v 2 ) chromoelectric dipole interaction and the lead- 
ing 0(v 2 ) spin-flip Fermi potential. Hence it is suppressed by the factor v A . Within the 
power counting assumed in this paper it only contributes in the N 4 LL approximation, far 
beyond the precision of our computation. Note that the nonperturbative contribution to 
the decay rates ratio is suppressed by a factor v 2 in comparison to the binding energy 
and decay rates, where the leading nonperturbative effect is due to chromoelectric dipole 
interaction. Thus the renormalization group improved result allows for very accurate 
theoretical evaluation of the spin ratio and, by using the available experimental data on 
the T meson as input, we can predict the production and annihilation rates of the yet 
undiscovered rjf, meson with high precision. In particular, we predict the rj^lS) decay 
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rate using the experimental value for the Y(IS') decay rate 

T(r] b (lS) -> 77) = 0.659 ± 0.089(th.)t^(<H) ± 0.015(exp.) keV , (23) 

where we have taken v = 1.295 GeV, the scale of minimal sensitivity, for the central 
value, the difference between the NLL and NNLL result for the theoretical error and 
a s (Mz) = 0.118 ± 0.003. The last error in Eq. ( !23l reflects the experimental error of 
r(T(15') — ► e + e") = 1.314 ±0.029 keV jl]. This value considerably exceeds the result for 
the absolute value of the decay width obtained in Ref. [B] on the basis of the full NLL 
analysis including the spin- independent part. This can be a signal of slow convergence of 
the logarithmic expansion for the spin-independent contribution which is more sensitive 
to the dynamics of the bound state and in particular to the ultrasoft contribution as it 
has been discussed above. On the other hand, the renormalon effects [45] could produce 
some systematic errors in the pure perturbative evaluations of the production/annihilation 
rates. The problem is expected to be more severe for the charmonium case discussed 
below. 

We would also like to remark that the one-loop result for v = m& overshoots the NNLL 
one by approximately 30%. This casts some doubts on the accuracy of the existing a s 
determination from the T(T — > light hadrons)/r(T — > e + e~) decay rates ratio, which 
gives a s (mb) = 0.177 ± 0.01, well below the "world average" value [4]. The theoretical 
uncertainty in the analysis is estimated through the scale dependence of the one-loop 
result. Our analysis of the photon mediated annihilation rates indicates that the actual 
magnitude of the higher order corrections is most likely quite beyond such an estimate 
and the theoretical uncertainty given in [4] should be increased by a factor of two. This 
brings the result for a s into la distance from the "world average" value. 

For the charmonium, the NNLO approximation becomes negative at an intermediate 
scale between a s m c and m c (c.f. Fig. [4]) and the use of the NRG is mandatory to get a 
sensible perturbative approximation. The NNLL approximation has good stability against 
the scale variation but the logarithmic expansion does not converge well. This is the 
main factor that limits the theoretical accuracy since the nonperturbative contribution 
is expected to be under control. For illustration, at the scale of minimal sensitivity, 
v = 0.645 GeV, we obtain 

Q^M^fl 1 (1 - 0.513 -0.340). (24) 

T(r] c (lS) -»■ 77) 3Q 2 C 

The central value of our NNLL result is la below the experimental value. The discrepancy 
may be explained by the large higher order contributions. This should not be surprising 
because of the rather large value of a s at the inverse Bohr radius of charmonium. For the 
charmonium hyperfine splitting, however, the logarithmic expansion converges well and 
the prediction of the NRG is in perfect agreement with the experimental data. Thus one 
can try to improve the convergence of the series for the production/annihilation rates by 
accurately taking into account the renormalon-related contributions. One point to note 
is that with a potential model evaluation of the wave function correction, the sign of the 
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NNLO term is reversed in the charmonium case [H]. At the same time the subtraction 
of the pole mass renormalon from the perturbative static potential makes explicit that 
the potential is steeper and closer to lattice and phenomenological potential models [46] . 
Therefore, the incorporation of higher order effects from the static potential may improve 
the agreement with experiment. In any case, if we estimate the theoretical uncertainty 
as the difference of the NNLL and the NLL result at the soft scale a s m c , the theoretical 
and experimental values agree within the error bars. 

6 Summary 

To conclude, we have resummed to all orders in the perturbative expansion the next- 
to-next-to-leading logarithms in the ratio of the one-photon mediated production and 
annihilation rate of the spin-triplet heavy quarkonium state to the two-photon one of the 
spin-singlet state. This constitutes the first complete NNLL result for the production and 
annihilation of heavy quarkonia. The use of the NRG improves the convergence of the 
perturbative series and stabilizes the result with respect to the scale variation. This allows 
us to obtain extremely accurate predictions for the tt threshold production. However, we 
still cannot draw a definite conclusion about the accuracy of the NNLL approximation to 
the spin independent part of the threshold production cross section, which could be worse 
due to the ultrasoft contribution. 

For the bottomonium case we found convergence of the logarithmic expansion at the 
physically motivated scale of inverse Bohr radius and very weak scale dependence of 
the NNLL approximation. Using these results we can give predictions for the yet to be 
discovered i]b meson two-photon production/annihilation rates. The magnitude of the 
higher order corrections questions the reported accuracy of a s determination from the 
T(T — > light hadrons)/r(T — > e + e~) ratio based on the one-loop theoretical analysis. 

For the charmonium annihilation the NRG is mandatory to bring physical sense to the 
perturbative result because the NNLO approximation becomes negative at an intermedi- 
ate scale. Though the NNLL approximation has very good stability with respect to the 
scale variation, the convergence of the logarithmic expansion is slow and essentially limits 
the accuracy of the perturbative prediction. Taking into account the large uncertainties 
of our result and the experimental data there is agreement and a lower value for the ratio 
seems to be favoured by our result. 

As a by-product of the heavy quarkonium analysis we have obtained novel results 
for positronium: the complete 0(a 4 ln 2 a) and 0(a 5 ln 3 a) corrections to the ratio of the 
ortho- and parapositronium decay rates. 
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Appendix A: Nonrelativistic renormalization group evo- 
lution of the spin-flip potential 

In order to obtain c s to NNLL accuracy one needs the matching coefficient of the spin-flip 

potential, -D52 s , to NLL order, which has been considered in Ref. [37] . In this Appendix 

we provide additional details. 

(2) 
The soft running of D s2 ' s to NLL approximation is discussed in detail in Ref. [37] . The 

NRG evolution of the spin-flip potential also includes the potential and ultrasoft running 
at NLL, which is described below. To compute the potential running, we inspect all op- 
erators that induce spin-dependent UV divergences in the time-independent perturbation 
theory contribution with one and two potential loops. They include 

(i) the tree- level 0(v A ) operators, 

(ii) the one-loop (D(a s v 3 ) operators, and 

(hi) the 0(v 2 , a s v) operators, 

which we in turn discuss in the following. In Appendix A. 4 we discuss the three-loop 
ultrasoft-potential contribution and, finally, in Appendix A. 5 we provide the NRG equa- 
tion for the potential running. 



A.l Tree-level potentials of order v A 

The operators relevant for our calculation read 

,(1)„(2) 



V. 



(-i) 



v« 



s 2 ,i 

(4) 



V t 



S 2 ,2 
(4) 



S 2 ,3 



V, 



(-1) 



S 2 A 



T^CpOl, 



—ttCfoi s 



c s' c s' 1 
Am 2 m 2 k 2 

r (l)J2) ■ 



cr l ■ (k x p)er 2 ■ (kx p) 



ycy \jr -p 



/2n2 



4mfm2 



kCf®. 



P -P 



,2 



2k' 



k z 

Am\m2 



<T\ ■ cr 2 — 



<j\ ■ k(T 2 ■ k 

~k 2 



>ix (p + p'))-((T 2 xk) + (1^2) 



tiCf&s 



J1)J2) 

u pp' u F 

•m\m 2 



(T! ■ (p + p') (a 2 ■ (p + p')k 2 



iP 



(25) 
(26) 

(27) 



p' 2 )cr 2 ■ k 



+ (!<-> 2) 



:28) 
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Here p and p' = p + k are the momentum of incoming and outgoing quark, we consider 
the general case with quark and antiquark of different masses nil and m 2 , and adopt the 
standard notations for the NRQCD coupling constants Cp, etc. (see, e.g. Ref. [47]). The 
operators in Eqs. ( l25l) - (l2~Tj) can be inferred from the analysis of the hyperffne splitting 
QED [25] with a trivial adjustment of the colour structure. Eq. (1251) corresponds to 
Eq. (13) of Ref. [25]. The potential (T2"6]) results from the expansion of the transverse 
gluon propagator in the energy transfer up to fc 2 with subsequent use of the Coulomb 
equation of motion [24] 

e (p'-^)\ (29) 

4mim 2 

It reproduces the retardation effect given by Eq. (32) of Ref. [23] . 

The potentials proportional to csCp appear both at tree and one-loop level. There are 
two cs NRQCD vertices that contribute to our calculation. One of them is proportional to 
AqA. It is responsible half of the contribution proportional to csCf in Eq. (130]) of the next 
subsection. The other involves the time derivative cfoA. To compute this contribution 
we perform a field redefinition in the NRQCD Lagrangian which in the lowest order is 
equivalent to using the equations of motion. In this way we obtain two new vertices. The 
first vertex is proportional to (p 2 /m)A and produces the potential in Eq. (f2TI) . which 
agrees with the first term in Eq. (22) of Ref. [25]. The second vertex is proportional 
to AqA and is responsible for the second half of the c,sCf contribution of the one-loop 
operator in Eq. (1311 . 

The potential (I2"gj) was not considered in Ref. [25] since in QED c p i p = 0(a) and 
thus V02 4 gives at most corrections of order a 3 to the hyperffne splitting. Moreover, in 
QED this potential contributes only to the NNLL running of the spin-flip operator as the 
coupling c p ip has no anomalous dimension. However, in QCD c p i p does run and we have 
to take it into account. 

There also exist 0(v A ) potentials including the product of cw 2 c f or cw^f- They are, 
however, proportional to (p • p')(cri • c 2 ) and (p 2 + p' )(<xi • cr 2 ), respectively, which do 
not produce potential divergences to the order of interest. 

To get the NLL approximation for the spin-flip potential we need the LL running of the 
0(v 4 ) potentials, which is determined by the soft running of the NRQCD coupling con- 
stants. They are either known [TB] or can be deduced using reparameterization invariance 



A. 2 One-loop potentials of order a s v 3 

The operators relevant for our calculation are 

, / (1) (2) (1) (2)\ 

V$ = -h™.) 2 C F (4C F -C A ) [S^Si^ + S^SiL \\k\<T l .cT 2 

24 I mi m 2 m i m 2 I 

1 r ( 1 ) r ( 2 ) 1 /„2, /2 k 2 \ 

H yKOis) Up^A ttt £& CT l " ^2 

24 mi 1^2 \k\ \ 2m r 2m r I 
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1 

96 1 



(1) (2)2 (1)2 (2) ' 



ircx x YC F C 



cv c 



F^A 



, °F u f Hi 



(30) 



where E is the two-particle energy. The Abelian piece reproduces Eq. (34) of Ref. [25] 
whereas the non- Abelian part is new. Note that the second line gives no NLL contribution 
as the term proportional to E develops no logarithm and the remaining two terms cancel 
each other. We need the NLL running of this potential, which is inherited through the 
LL soft running of the NRQCD coupling constants. 



A. 3 Double insertion of 0(v 2 ,a s v) operators 

The basis of the 0(v 2 , a s v) operators is well known from the lower order calculations and 
can, e.g. be found in Refs. [48,27] . To get the NLL contribution to the spin-flip potential 
from the double insertion of these operators we need the LL soft and ultrasoft running 
of the 0(v 2 ) terms and NLL soft and ultrasoft running of 0(a s v) operators, which are 
known for the equal mass case [16]. For the non-equal mass case, most of the results can 
be trivially generalized (c.f. Ref. [16]) except for D ds , which depends on the four-fermion 
couplings of NRQCD. Among these couplings only for d vs a nontrivial change is necessary. 
The corresponding NRG equation can be found in Ref. [16] and its solution reads 



dvaiy) = d vs {v h ) - ot s (v h ) 



x- 



-ot a {v h ) 

A) 



\ Cf - 3^ - \c A (™ + ?*)] ^ (; 

2 4 A Vm 2 mj\ A V 

/mi m 2 \ 27C\ vr J 5C A + AT F m 

\^~ 2 + ~^J 9C A + 8T F m 2/3 \4C A + 4T Fni 



A) 



- 1 



A) - 2C A 



z 



O -2C A 



"I + 



+ 



-7(7/i + 32C A C F 



C A + 1QC F - 8T F rii 
2(C A - 2T Fni ) 
4C A T F ni + 32C F T F ni 



(> - l) 



x- 



4(C A + T F nO(2T F ^-C A ) 
3A, 



3/3 + AT F m - 2C A 



/3 +4T F n,/3-2C A /3 



+ 



8T F m 



9C> 



0o 



X 



y 3o 



A) - 2C, 



-1 



z 



P -2C A 



1 



1 + 



6Ao 



6Ao - 13C A 



y 0o- 



nc A 



20 
13" 

-1 



32 Cf 
13^ 



(31) 



A. 4 Three-loop ultrasoft-potential running 

Another NLL contribution to the spin-flip potential comes from the three-loop NLL 
ultrasoft-potential running of a 0(v 2 ) spin-flip operator. The contribution is related to 
the chromomagnetic dipole gluon exchange which generate the potential of the following 
form in the coordinate representation 



V, 



(2) 



5 2 ,l/r3 



4CVS 2 
3mim2 



V (V - V S ) 2 D 



(2) 
5 2 ,l/r 3 ' 



(32) 
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where 



v n 



<^F] — r ? 
\r\ 




T" C ') 


\r 



(33) 



The ultrasoft running of the potential is determined by the NRG equation 



dV, 



(2) 



5 2 ,l/r 3 



V„ 



d In z/„ 



C-FK^us) 



a s (v us ) 2C F 

7T Zm\V(l2 



Vo{Vo - v s f 



and the solution for the Wilson coefficient reads 



D (2) (v )- — 



' g s {vh) 

a s (v us ) 



2C A /0o 



« S (^) 
aJl/\r\ 



2C A /0o' 



(34) 



(35) 



The potential ( 1321) is singular at \r\ — > and should be understood as the Fourier trans- 
form of \\i(\k\/v p ). 



A. 5 NRG equation for potential running 

With all the relevant operators at hand it is straightforward to derive the NRG equation 
for the NLL potential running of the spin-flip potential. It reads 



dD { sl 
dlnu 



-2m: 



rri'l m^J ' m\m2 



(2) n (2) 



+- (n {2) 



w (2) D (2) _ _5_ / D (2) n2 



C F a Vs \2D%D% 

c a c f d^d£1 



+ U C F m r f C ^ C F + CF C M ) - ^ICl^-cfcf 



ml 



rrin 



771x1712 

& „(1) 



1z / \^} K^J 



— afC 



1 ' " F ^^ C5 ° 5 asCpmr \ Cf ml ril F 



(2) (1) 

+$a£C F {4C F - C A )m r I cT— + —cf 

2 \ 7712 77l\ 



+ -a 3 C F CAm r c 



(1)C F 



.(2) 2 (1) 2 



+ 



F ' C F 

m 2 mi 



^\ ( 2)\ C\{C A -2C F )a* ( 2 ) 



D 



5 2 ,l/r3 ' 



(36) 



where m r = mim2/(mi+m2) is the reduced mass. The notation for the Wilson coefficients 
is adopted from Ref. [16], where also explicit results can be found. The first two lines 
of Eq. (136]) correspond to the double insertion of the 0(v 2 , a s v) operators discussed in 
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Appendix A. 3. Note that the contribution proportional to [Dg* J disagrees with the 
result of Ref. [H]. To reproduce the 0(a 2 ) corrections to the positronium hyperfine 
splitting [oT)ll2o] in Ref. |IH] the factor 9 of the corresponding Wilson coefficient U 2 should 
be changed to 5. This was corrected in Ref. [5T]. The expressions in the third and forth 
line of Eq. (136|) follow from the results for Appendix A.l. The fifth and the first term 
of the sixth line are obtained from Appendix A. 2. Finally, the last term in Eq. (1361) 
corresponds to the ultrasoft-potential running discussed in Appendix A. 4. 

Appendix B: Analytical result for the NNLL contri- 
bution to c s 

The analytical results for A{ and fi(z) of Eq. ( fl9l) read 

1C A z&\ 



2 



2p Q -(25C A )/e 



/l( , )= ^o-2^ 3F2 ^3_^ )3 _2CU ;4 _2CU )4 _2^ 
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f 8 (z) = /0-(13^)/6 j fg{z) = z ?o-2C A ^ /loW = /0 ; f n{z)=Hz)} 

f 12 (z) = \n 2 (z) , f 13 (z) = In (2 - z*>) , / M (z) = ^" 2 ^ In (2 - z*>) , 
/ 15 (*)=/°ln(2-A) , 

/„(z) = 2 *W/V*-*^2 f-^, 2 - ^, 2 - ^; 3 - ^, 3 - ^; ^) , 

\ Po Po Po Po Po l j 

f 17 (z) = l, f 18 ( z )=z 2 ^, f w (z) = IA 2 (f>/2). (37) 
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where B z is the incomplete beta-function, and p F q is the hypergeometric function. 
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Appendix C: 0(a 4 ln 2 a) and 0(oP\y? a) corrections to 
the ratio of ortho- and parapositronium decay rates 

We can also obtain novel results for the ratio of the ortho- (oPs) and parapositronium 
(pPs) decay rates by considering the Abelian limit of the previous computation and 
taking into account the modifications of the NRQCD matching coefficients due to the 
annihilation terms. The series for the decay rates ratio reads 



TjoPs -> fry) 4(tt 2 - 9) 

T{pPs -» 27) ~ 9tt 
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36 07T J 

where A = 10.286606(10), B = 44.87(26) and B v = 5.1243(33). The first four correction 
terms in the curly brackets can be extracted from know results (see e.g. [52]). The last 
two terms are our new results. 
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